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ABSTRACT 

We review uses of the generalized-ensemble algorithms for free-energy calculations in 
protein folding. Two of the well-known methods are multicanonical algorithm and replica- 
exchange method; the latter is also referred to as parallel tempering. We present a new 
generalized-ensemble algorithm that combines the merits of the two methods; it is referred 
to as the replica-exchange multicanonical algorithm. We also give a multidimensional 
extension of the replica-exchange method. Its realization as an umbrella sampling method, 
which we refer to as the replica-exchange umbrella sampling, is a powerful algorithm that 
can give free energy in wide reaction coordinate space. 

1 Introduction 

Over the past three decades, a number of powerful simulation algorithms have been intro- 
duced to the protein folding problem (for reviews see, e.g., Refs. For many years, 
the emphasis has been placed on how to find the global-minimum-energy conformation 
among a huge number of local-minimum states. For complete understanding of protein 
folding mechanism, however, the global knowledge of the configurational space is required, 
including the intermediate and denatured states of proteins. For this purpose, free-energy 
calculations are essential. 

We have been advocating the uses of generalized-ensemble algorithms as the methods 
that meet the above requirements (for reviews see, e.g., Refs. [§, ||). In this method 
each state is weighted by a non-Boltzmann probability weight factor so that a random 
walk in potential energy space may be realized. The random walk allows the simulation 
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to escape from any energy barrier and to sample much wider configurational space than 
by conventional methods. Monitoring the energy in a single simulation run, one can 
obtain not only the global-minimum-energy state but also canonical ensemble averages 
as functions of temperature by the single-histogram || and/or multiple-histogram J/| |^| 
reweighting techniques (an extension of the multiple-histogram method is also referred to 
as weighted histogram analysis method (WHAM) ||). 

Three of the most well-known generalized-ensemble methods are perhaps multicanon- 
ical algorithm (MUCA) || [Kj, simulated tempering (ST) [|ll|, [12|], and replica- exchange 
method (REM) ]T3|, [TJj]. (MUCA is also referred to as entropic sampling JT5I , [T(| and 
adaptive umbrella sampling |l7j . ST is also referred to as the method of expanded ensem- 
ble W% - REM is also referred to as parallel tempering [18|. Details of literature about 



REM and related algorithms can be found in a recent review |19| . ) Since MUCA was first 
introduced to protein folding problem f2(| , various generalized-ensemble algorithms have 
been used in many applications in protein and related systems (see Ref. || and references 
therein). In particular, free-energy calculations in protein folding by generalized-ensemble 
algorithms were explored in Refs. |2~T|, E2J. 



REM has been drawing much attention recently because the probability weight factors 
are essentially known a priori, whereas they are not in most of other generalized-ensemble 
algorithms (and have to be determined by a tedius procedure). In REM a number of 
non-interacting copies (or replicas) of the original system at different temperatures are 
simulated independently and simultaneously by the conventional Monte Carlo (MC) or 
molecular dynamics (MD) method. Every few steps, pairs of replicas are exchanged with 
a specified transition probability. 

We have worked out the details for the replica-exchange molecular dynamics algorithm 
p3| (see also Ref. ||24[ ). We have also developed a multidimensional replica- exchange 
method (MREM) JZjJ (see also Refs. g|, HI)- In M R EM we showed that REM is not 
limited to tempering (or temperature exchange) and that we can also exchange parameters 
in the potential energy. Namely, pairs of replicas with different temperatures and/or dif- 
ferent parameters of the potential energy are exchanged during the simulation. Important 
applications of MREM are free-energy calculations. 

The umbrella sampling method |28] and free energy perturbation method, which is a 
special case of umbrella sampling, have been widely used to calculate the free energies 
in chemical processes |28| - [37|. Although the effectiveness of the umbrella sampling 
method is well known, its successful implementation requires a careful fine tuning. Various 
generalizations of the umbrella sampling method have thus been introduced to sample the 
potential energy surface more effectively. The X-dynamics - |HD[ is such an example, 



where the coupling parameter A is treated as a dynamical variable. Another example is the 



multicanonical WHAM which combines the umbrella sampling with multicanonical 
algorithm. We have developed yet another generalization of the umbrella sampling method 
(we refer to this method as replica- exchange umbrella sampling (REUS)), which is based 
on the multidimensional extension of the replica-exchange method [ p5[| . 

REM is very effective and has already been used in many applications in protein 
systems (see Ref. []5| and references therein). However, REM also has a computational 
difficulty: As the number of degrees of freedom of the system increases, the required num- 
ber of replicas also greatly increases, whereas only a single replica is simulated in MUCA 
or ST. This demands a lot of computer power for complex systems. Our solution to this 
problem is: Use REM for the weight factor determinations of MUCA or ST, which is much 



2 



simpler than previous iterative methods of weight determinations, and then perform a long 
MUCA or ST production run. The first example is the replica- exchange multicanonical al- 
gorithm (REMUCA) In REMUCA, a short replica-exchange simulation is performed, 
and the multicanonical weight factor is determined by WHAM [0, ||. Another example of 
such a combination is the replica- exchange simulated tempering (REST) |[43| . In REST, 
a short replica-exchange simulation is performed, and the simulated tempering weight 
factor is determined by WHAM [|7|, |J. We have introduced a further extension of RE- 
MUCA, which we refer to as multicanonical replica-exchange method (MUCAREM) j42fl . 
In MUCAREM, the multicanonical weight factor is first determined as in REMUCA, and 
then a replica-exchange multicanonical production simulation is performed with a small 
number of replicas. 

In this article, we first describe the multidimensional replica-exchange method, a par- 



ticular realization of which is the replica-exchange umbrella sampling ||25|| . We then 



present the replica-exchange multicanonical algorithm |]l2"fl . The effectiveness of these 
methods is tested with short peptide systems. 



2 Methods 

2.1 Multidimensional Replica- Exchange Method 



We first review the original replica- exchange method (REM) | | (see Ref. |23| for 
details of the molecular dynamics version). 

We consider a system of N atoms with their coordinate vectors and momentum vectors 
denoted by q = {q\, ■ ' ' iQn} an d V = {Pi> • • • ,Pn}j respectively. The Hamiltonian 
H(q,p) of the system is the sum of the kinetic energy K(p) and the potential energy 
E{q): 

H{q,p)=K{p) + E{q) . (1) 

In the canonical ensemble at temperature T each state x = (q,p) with the Hamiltonian 
H(q,p) is weighted by the Boltzmann factor: 

W T (x) = e ~ pH ^ , (2) 

where the inverse temperature (3 is defined by /3 = 1/ksT (fc# is the Boltzmann constant). 

The generalized ensemble for REM consists of M non-interacting copies (or, repli- 
cas) of the original system in the canonical ensemble at M different temperatures T m 
(m = 1,---,M). We arrange the replicas so that there is always exactly one replica 
at each temperature. Then there is a one-to-one correspondence between replicas and 
temperatures; the label i (i = 1, • ■ ■ , M) for replicas is a permutation of the label m 
(m = 1, • • • , M) for temperatures, and vice versa: 

= i(m) = f(m) , , . 

= m(t) = f-\i) , [6) 

where f(m) is a permutation function of m and / _1 («) is its inverse. 

Let X = (xi^\ ■ ■ ■ , £jv/ M ^} = {^m(i), • • ■ , ^m(M) } stand for a "state" in this gener- 
alized ensemble. Here, the superscript and the subscript in ajM label the replica and the 
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temperature, respectively. The state X is specified by the M sets of coordinates and 
momenta p^ of X atoms in replica i at temperature T m : 



(4) 



Because the replicas are non-interacting, the weight factor for the state X is given by the 
product of Boltzmann factors for each replica (or at each temperature): 



WW*) =exp|-^/3 m(i) ^(gH,pH)| , 

( M 

= exp I - ]T /3 m i? (glWl./Hl 



(5) 



m=l 



where i(m) and m(i) are the permutation functions in Eq. (^). 

We now consider exchanging a pair of replicas. Suppose we exchange replicas % and j 
which are at temperatures T m and T n , respectively: 



The exchange of replicas can be written in more detail as 



',-} 



x% = g tJ ,p 



(6) 



(7) 



where the momenta are uniformly rescaled according to [23 



T 



m 



V" =\\—V 
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In order for this exchange process to converge towards the equilibrium distribution 
based on Eq. (|5J), it is sufficient to impose the detailed balance condition on the transition 
probability w(X — > X'): 



W REM (X) w{X - X') = W REM (X>) w(X> - X) 

From Eqs. ©, ©, ®, and @, we have 

iu(X -> X') 



w(X' -> X) 



exp (—A) 



where 



= (/3 m - /? n ) (E (g^) - E (V* 
This can be satisfied, for instance, by the usual Metropolis criterion: 



w{X -*')=«, (xg 



[j] 



1 , for A < 

exp (-A) , for A > 



(9) 
(10) 

(11) 
(12) 

(13) 



Note that because of the velocity rescaling of Eq. (§) the kinetic energy terms are cancelled 
out in Eqs. (O) (and (|T2"D) and that the same criterion, Eqs. ([II]) and (|TB|), which was 
originally derived for Monte Carlo algorithm |13], [L4| is recovered 



A simulation of the replica-exchange method (REM) |T3|, |Tj] is then realized by alter- 
nately performing the following two steps: 

1. Each replica in canonical ensemble of the fixed temperature is simulated simultaneously 
and independently for a certain MC or MD steps. 



2. A pair of replicas, say and x$, are exchanged with the probability w {x\^ 
in Eq. (0). 



r b1 



In the present work, we employ molecular dynamics algorithm for Step 1. Note that 
in Step 2 we exchange only pairs of replicas corresponding to neighboring temperatures, 
because the acceptance ratio of the exchange decreases exponentially with the difference 
of the two /3's (see Eqs. ([12]) and (|T3|)). Note also that whenever a replica exchange is 
accepted in Step 2, the permutation functions in Eq. (|3|) are updated. 

The method is particularly suitable for parallel computers. Because one can minimize 
the amount of information exchanged among nodes, it is best to assign each replica to each 
node (exchanging pairs of temperature values among nodes is much faster than exchanging 
coordinates and momenta). This means that we keep track of the permutation function 
m(i; t) = t) in Eq. (|[) as a function of MD step t throughout the simulation. 

The major advantage of REM over other generalized-ensemble methods such as mul- 



ticanonical algorithm [|| [L0| and simulated tempering [|TT|, [12|] lies in the fact that the 
weight factor is a priori known (see Eq. (||) ) , whereas in the latter algorithms the deter- 
mination of the weight factors can be very tedius and time-consuming. A random walk 
in "temperature space" is realized for each replica, which in turn induces a random walk 
in potential energy space. This alleviates the problem of getting trapped in states of 
energy local minima. In REM, however, the number of required replicas increases as the 
system size N increases (according to \/N) This demands a lot of computer power 
for complex systems. 

We now present our multidimensional extension of REM, which we refer to as mul- 
tidimensional replica- exchange method (MREM). The crucial observation that led to the 
new algorithm is: As long as we have M non-interacting replicas of the original system, 
the Hamiltonian H(q,p) of the system does not have to be identical among the replicas 
and it can depend on a parameter with different parameter values for different replicas. 
Namely, we can write the Hamiltonian for the 2-th replica at temperature T m as 

H m ^\p^) = K(p®) + E Xm (q®) , (14) 

where the potential energy E\ m depends on a parameter A m and can be written as 

E Xm (qtt) = E ( q M) + \ m V(qM) . (15) 

This expression for the potential energy is often used in simulations. For instance, in 



umbrella sampling p8 |, Eo(q) and V(q) can be respectively taken as the original potential 
energy and the "biasing" potential energy with the coupling parameter A m . In simulations 
of spin systems, on the other hand, E (q) and V(q) (here, q stands for spins) can be 
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respectively considered as the zero-field term and the magnetization term coupled with 
the external field X m . 

While replica i and temperature T m are in one-to-one correspondence in the original 
REM, replica i and "parameter set" A m = (T m , X m ) are in one-to-one correspondence 
in the new algorithm. Hence, the present algorithm can be considered as a multidimen- 
sional extension of the original replica-exchange method where the "parameter space" is 
one-dimensional (i.e., A m = T m ). Because the replicas are non-interacting, the weight 
factor for the state X in this new generalized ensemble is again given by the product of 
Boltzmann factors for each replica (see Eq. (H)): 



Wmrem(AT) = exp | -J2 Pm(i)H m (i) (V i] , P W ) | 

( M 



(16) 



m=l 



where i(m) and m(i) are the permutation functions in Eq. (0). Then the same derivation 
that led to the original replica-exchange criterion follows, and the transition probability 
of replica exchange is given by Eq. (|T3|), where we now have (see Eq. (O)) 



A 



An [E\ m ( 



[.)} 



Pn (E Xn (q 



,b1 



(17) 



Here, E\ m and E Xn are the total potential energies (see Eq. (|I5|) ). Note that we need to 
newly evaluate the potential energy for exchanged coordinates, E\ m (q^') and E\ n (q^), 
because E\ m and E\ n are in general different functions. 

For obtaining the canonical distributions, the weighted histogram analysis method 
(WHAM) m is particularly suitable. Suppose we have made a single run of the present 
replica-exchange simulation with M replicas that correspond to M different parameter 
sets A m = (T m , A m ) (m = 1, ■ ■ • , M). Let N m (E Q , V) and n m be respectively the potential- 
energy histogram and the total number of samples obtained for the m-th parameter set 
A m . The WHAM equations that yield the canonical probability distribution Pt^Eq, V) 
with any potential-energy parameter value A at any temperature T = l/ks/3 are then 
given by ML § 



Pt,x(E ,V) 



M 



J2 g- 1 N m (E ,V) 



m=l 



J2 9r, 



n, 



m=l 



and 



3 /m ftmE\ r 



,(E ,V) 



-PE X 



(19) 



e -/« = £ p Tn 

E ,V 

Here, g m = 1 + 2r m , and r m is the integrated autocorrelation time at temperature T m with 
the parameter value \ m . Note that the unnormalized probability distribution Pt,\(Eq, V) 
and the "dimensionless" Helmholtz free energy f m in Eqs. (18) and (19) are solved self- 
consistently by iteration [J71 

We can use this new replica-exchange method for free energy calculations. We first 
describe the free-energy perturbation case. The potential energy is given by 



E\(q) — E I (q) + X (E F (q) — E I (q)) 



(20) 
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where Ej and Ep are the potential energy for a "wild-type" molecule and a "mutated" 
molecule, respectively Note that this equation has the same form as Eq. (|i"5|). 

Our replica-exchange simulation is performed for M replicas with M different values 
of the parameters A m = (T m ,A m ). Since E\ =0 (q) = Ej(q) and E\=i(q) = E F (q), we 
should choose enough A m values distributed in the range between and 1 so that we 
may have sufficient replica exchanges. From the simulation, M histograms N m (Ei, Ep — 
Ej), or equivalently N m (Ej, Ep), are obtained. The Helmholtz free energy difference of 
"mutation" at temperature T, AE = F\ = i — F\ =Q , can then be calculated from 



exp(-pAF) 



J2 Pt,x=i(E!,E f ) 

*T,\=1 _ Ej,Ep 

?t,a=o ~ £ Pt,x=o(E t , Ep) 

Ej,E F 



[21] 



where Pt,\(Ei, Ep) are obtained from the WHAM equations of Eqs. (jl8D and (J19|) . 

We now describe another free energy calculations based on MREM applied to umbrella 
sampling pgfl , which we refer to as replica- exchange umbrella sampling (REUS). The 
potential energy is a generalization of Eq. (|T5D and is given by 



E x (q)=E (q)+J2^ ) V e (q) 



(22) 



i=i 



where E (q) is the original unbiased potential, Vi(q) (£ = 1, • • • , L) are the biasing (um- 
brella) potentials, and are the corresponding coupling constants (A = (A^, . . . ; A^)). 
Introducing a "reaction coordinate" £, the umbrella potentials are usually written as har- 
monic restraints: 



Ve(q) = h Uq) - d t f , {£ 



(23) 



where dg are the midpoints and kg are the strengths of the restraining potentials. We 
prepare M replicas with M different values of the parameters A m = (T m ,A m ), and the 
replica-exchange simulation is performed. Since the umbrella potentials Ve(q) in Eq. (f2"3"| ) 
are all functions of the reaction coordinate £ only, we can take the histogram N m (E ,^) 
instead of N m (E , V x , ■ ■ ■ , V L ). The WHAM equations of Eqs. (0) and (TTJ) can then be 
written as 



M 



^g- 1 N m (E ,Z) 



m=l 



M 



^ y 9m 



II, 



frn—ftmE \ 
3 An 



m=l 



and 



(24) 



(25) 



The expectation value of a physical quantity A is now given by 



< A > 



Y i A(E ,OP TX (E ,0 



(26) 
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The potential of mean force (PMF), or free energy as a function of the reaction coor- 
dinate, of the original, unbiased system at temperature T is given by 



where {0} = (0, •••,0). 



E 



(27) 



2.2 Replica- Exchange Multicanonical Algorithm 



We first briefly review the multicanonical algorithm || |10[- Because the coordinates q 
and momenta p are decoupled in Eq. ([!]), we can suppress the kinetic energy part and can 
write the Boltzmann factor as 

W T (x) = W T (E) = e~ pE . (28) 

The canonical probability distribution of potential energy Pt{E) is then given by the 
product of the density of states n(E) and the Boltzmann weight factor Wt{E): 

Pt(E) oc n(E)W T (E) . (29) 

In the multicanonical ensemble (MUCA) || [T(J, on the other hand, each state is 
weighted by a non-Boltzmann weight factor W mn (E) (which we refer to as the multi- 
canonical weight factor) so that a uniform energy distribution P mn (E) is obtained: 

P nm (E) oc n(E)W mn {E) = constant . (30) 

The flat distribution implies that a free random walk in the potential energy space is real- 
ized in this ensemble. This allows the simulation to escape from any local minimum-energy 
states and to sample the configurational space much more widely than the conventional 
canonical MC or MD methods. 

From the definition in Eq. ( |30| ) the multicanonical weight factor is inversely propor- 
tional to the density of states, and we can write it as follows: 

W mu {E) = e -/WW*;*b) = J_ (31) 

n(E) 

where we have chosen an arbitrary reference temperature, T = 1//cbA), an d the "multi- 
canonical potential energy" is defined by 

E mu (E;T ) = k B T lnn(E) = T S(E) . (32) 

Here, S(E) is the entropy in the microcanonical ensemble. 

A multicanonical Monte Carlo simulation is performed, for instance, with the usual 
Metropolis criterion: The transition probability of state x with potential energy E to 
state x' with potential energy E' is given by 

1 ] \ exp (-(3 AE mu ) , for AE mu > , [66) 
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where 

Ai? mu = E mu (E'; T ) — E mu (E; T ) . (34) 

The molecular dynamics algorithm in multicanonical ensemble also naturally follows from 
Eq. fl3~T|), in which the regular constant temperature molecular dynamics simulation (with 
T = To) is performed by solving the following modified Newton equation: [44, (i^, [T7|| 



dE mu (E; T ) dE mu (E; T ) 
Pk ~ t\ k ~ BE Tk ' [ } 



where f k is the usual force acting on the fc-th atom (k = 1, • ■ ■ , N). From Eq. (p2|) this 
equation can be rewritten as 

T 

Pk = fk , (36) 

where the following thermodynamic relation gives the definition of the "effective temper- 
ature" T(E): 

aS < E » - (37) 



dE 
with 



E=E a T ( E ' 



E a = <E> T{Ea) . (38) 

The multicanonical weight factor is usually determined by iterations of short trial 
simulations. The details of this process are described, for instance, in Refs. jit], [46| . 
However, the iterative process can be non-trivial and very tedius for complex systems. 

After the optimal multicanonical weight factor is determined, one performs a long 
multicanonical simulation once. By monitoring the potential energy throughout the sim- 
ulation, one can find the global-minimum-energy state. Moreover, by using the obtained 
histogram N mn (E) of the potential energy distribution P mu (E), the expectation value of 
a physical quantity A at any temperature T = l/ksP can be calculated from 



< A >t = E ^ — _ pE , (39) 



E 



where the best estimate of the density of states n(E) is given by the single-histogram 
reweighting techniques (see Eq. (|30|)) ||: 

- tS! • (4o) 

The replica- exchange multicanonical algorithm (REMUCA) overcomes both the 
difficulties of MUCA (the multicanonical weight factor determination is non-trivial) and 
REM (a lot of replicas, or computation time, is required). In REMUCA we first perform 
a short REM simulation (with M replicas) to determine the multicanonical weight factor 
and then perform with this weight factor a regular multicanonical simulation with high 
statistics. The first step is accomplished by the weighted histogram analysis method 0. Q . 
Let N m (E) and n m be respectively the potential-energy histogram and the total number 
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of samples obtained at temperature T m = 1/ &b An of the REM run. The density of states 
n(E) is then given by solving the following WHAM equations [j], |^|: 



M 

£ 97n N m (E) 
<E) = , (41) 

y m 'I'm c 

m=l 

where 

e -/m = n{E) e~ 0mE . (42) 

E 

Once the estimate of the density of states is obtained, the multicanonical weight factor 
can be directly determined from Eq. ([31]) (see also Eq. (|32|) ). Actually, the multicanonical 
potential energy E mu (E;T ) thus determined is only reliable in the following range: 

E\ < E < Em , (43) 

where 

| E x = <E> Tl , 
\ E M = < E > Tm , 

and Ti and Ta/ are respectively the lowest and the highest temperatures used in the REM 
run. Outside this range we extrapolate the multicanonical potential energy linearly: 



dE mu (E; T 




{E - E) + E mn (E i; T ) , for E < E 1 



ej£(E) = { E mu (E;T ) , E El for E x < E < E M , (45) 

(E - E M ) + E mu {E M ] T ) , for E > E M • 

A long multicanonical MD run is then performed by solving the Newton equations in 
Eq. (|35D into which we substitute £^(E) of Eq. (|^). Finally, the results are analyzed 
by the single-histogram reweighting techniques as described in Eq. fl4"0D (and Eq. (PP|)). 
We remark that our multicanonical MD simulation here actually results in a canonical 
simulation at T = T\ for E < E\, a multicanonical simulation for E\ < E < Em, and a 
canonical simulation at T = Tm for E > Em (a detailed discussion on this point is given 
in Ref. [§]). Note also that the above arguments are independent of the value of T , and 
we will get the same results, regardless of its value. 

Since the WHAM equations are based on histograms, the density of states n(E), or 
the multicanonical potential energy E mu (E]T ), will be given in discrete values of the 
potential energy E. For multicanonical MD simulations, however, we need the derivative 
of E mn (E;T ) with respect to E (see Eq. (^5])). We thus introduce some smooth function 
to fit the data. It is best to fit the derivative d gm ^' T °) directly rather than E mu (E;T ) 
itself. For this we recall the Newton equation of Eq. ( |3T)D and the thermodynamic relation 
of Eqs. ([37]) and ([38]). The effective temperature T(E), or the derivative d£muCE,T ) ^ can 
be obtained by fitting the inverse of Eq. (BH) by a smooth function, where < E >t(e) is 



calculated from Eq. ( p9[) by solving the WHAM equations of Eqs. fl41]) and (0). Given 
its derivative, the multicanonical potential energy can be obtained by integration: 

To) = To k -dE~ dE = To k TjE) • (46) 
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We remark that the same equations were used to obtain the multicanonical weight factor 
in Ref. |4"7|] , where < E >t(e) was estimated by simulated annealing instead of REM. 

Finally, although we did not find any difficulty in the case of protein systems that we 
studied, a single REM run in general may not be able to give an accurate estimate of the 
density of states (like in the case of a strong first-order phase transition |13[). In such a case 
we can still greatly simplify the process of the multicanonical weight factor determination 
by combining the present method with the previous iterative methods |PH [16f . 

The formulation of REMUCA is simple and straightforward, but the numerical im- 
provement is great, because the weight factor determination for MUCA becomes very 
difficult by the usual iterative processes for complex systems. 



3 RESULTS 

We now present some examples of the simulation results by the algorithms described in 
the previous section. Short peptide systems were considered. 

For molecular dynamics simulations, the force-field parameters were taken from the 
all-atom versions of AMBER ED[. The computer code developed in Refs. |5"0|, 



which is based on PRESTO |plj| , was used. The unit time step was set to 0.5 fs. The 
temperature during the canonical MD simulations was controlled by the constraint method 
1 52], [53|]. Besides gas phase simulations, we have also performed MD simulations with a 



distance-dependent dielectric, e = r, and with explicit water molecules of TIP3P model 
|54[. 

As described in detail in the previous section, in generalized-ensemble simulations 
and subsequent analyses of the data, potential energy distributions have to be taken as 
histograms. For the bin size of these histograms, we used the values ranging from 0.5 to 
2 kcal/mol, depending on the system studied. 

The first example is a penta peptide, Met-enkephalin, whose amino-acid sequence 
is: Tyr-Gly-Gly-Phe-Met. This peptide in gas phase was studied with the force field 



of AMBER in Ref. [48| by the replica-exchange MD simulation |[23|| . We made an MD 
simulation of 2 x 10 6 time steps (or, 1.0 ns) for each replica, starting from an extended 
conformation. We used the following eight temperatures: 700, 585, 489, 409, 342, 286, 
239, and 200 K, which are distributed exponentially, following the annealing schedule 



of simulated annealing simulations p5| . As is shown below, this choice already gave an 



optimal temperature distribution. The replica exchange was tried every 10 fs, and the 
data were stored just before the replica exchange for later analyses. 

As for expectation values of physical quantities at various temperatures, we used 



the weighted histogram analysis method of Eqs. fl4T|) and (42). We remark that for 
biomolecular systems the integrated autocorrelation times r m in the reweighting formulae 
(see Eq. (|4lD) can safely be set to be a constant |§, and we do so throughout the analyses 
in this section. 

In Figure 1 the time series of temperature exchange (a) and the total potential energy 
(b) for one of the replicas are shown. We do observe random walks in both temperature 
space and potential energy space. Note that there is a strong correlation between the 
behaviors in Figures 1(a) and 1(b). 

In Figure 2 the canonical probability distributions obtained at the chosen eight tem- 
peratures from the replica-exchange simulation are shown. We see that there are enough 
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Figure 1: Time series of (a) temperature exchange and (b) the total potential energy 
for one of the replicas from a replica-exchange MD simulation of Met-enkephalin in gas 
phase. 
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overlaps between all pairs of distributions, indicating that there will be sufficient numbers 
of replica exchanges between pairs of replicas. 




Figure 2: The canonical probability distributions of the total potential energy of Met- 
enkephalin in gas phase obtained from the replica-exchange MD simulation at the eight 
temperatures. The distributions correspond to the following temperatures (from left to 
right): 200, 239, 286, 342, 409, 489, 585, and 700 K. 



We further compare the results of the replica-exchange simulation with those of a 
single canonical MD simulation (of 1 ns) at the corresponding temperatures. In Figure 
3 we compare the distributions of a pair of main-chain dihedral angles ((f), ip) of Gly-2 
at two extreme temperatures (T = 200 K and 700 K). While the results at T = 200 K 
from the regular canonical simulation are localized with only one dominant peak, those 
from the replica-exchange simulation have several peaks (compare Figures 3(a) and 3(b)). 
Hence, the replica-exchange run samples much broader configurational space than the 
conventional canonical run at low temperatures. The results at T = 700 K (Figures 3(c) 
and 3(d)), on the other hand, are similar, implying that a regular canonical simulation 
can give accurate thermodynamic quantities at high temperatures. 

In Figure 4 we show the average total potential energy as a function of temperature. 
As expected from the results of Figure 3, we observe that the canonical simulations at low 
temperatures got trapped in states of energy local minima, resulting in the discrepancies 
in average values between the results from the canonical simulations and those from the 
replica-exchange simulation. 

We now present the results of replica- exchange umbrella sampling (REUS) [25|| . The 
system of a blocked peptide, alanine-trimer, was studied. Since the thermodynamic be- 
havior of this peptide was extensively studied by the conventional umbrella sampling |JT 



it is a good test case to examine the effectiveness of the new method. The force field pa- 
rameters were taken from the all-atom version of AMBER [4~8"| with a distance-dependent 
dielectric, e = r, which mimics the presence of solvent. We made an MD simulation of 
4 x 10 6 time steps (or, 2.0 ns) for each replica, starting from an extended conformation. 
The data were stored every 20 steps (or, 10 fs) for a total of 2 x 10 5 snapshots. 

In Table 1 we summarize the parameters characterizing the replicas for the simula- 
tions performed. They are one original replica-exchange simulation (REM1), two replica- 
exchange umbrella sampling simulations (REUS1 and REUS2), and two conventional 
umbrella sampling simulations (US1 and US2). 
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Figure 3: Distributions of a pair of main-chain dihedral angles (</>, ip) of Gly-2 for: (a) 
T = 200 K from a regular canonical MD simulation, (b) T = 200 K from the replica- 
exchange MD simulation, (c) T = 700 K from a regular canonical MD simulation, and 
(d) T = 700 K from the replica-exchange MD simulation. 




200 300 400 500 600 700 
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Figure 4: Average total potential energy of Met-enkephalin in gas phase as a function of 
temperature. The solid curve is the result from the replica-exchange MD simulation and 
the dots are those of regular canonical MD simulations. 
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Table 1: Summary of Parameters in US, REM, and REUS Simulations 
Run a M 5 iVr 5 Temperature [K] T b d t [A] (k e [kcal/mol-A 5 



2l\c 



REM1 16 16 200, 229, 262, 299, 

342, 391, 448, 512, 
586, 670, 766, 876, 
1002, 1147, 1311, 
1500 

REUS1, US1 14 1 300 14 0.0 (0.0) d , 1.8 (1.2), 2.8 (1.2), 3.8 (1.2), 

4.8 (1.2), 5.8 (1.2), 6.8 (1.2), 7.8 (1.2), 
8.8 (1.2), 9.8 (1.2), 10.8 (1.2), 11.8 (1.2), 
12.8 (1.2), 13.8 (1.2) 

REUS2, US2 16 4 250, 315, 397, 500 4 0.0 (0.0), 7.8 (0.3), 10.8 (0.3), 13.8 (0.3) 



a REM, REUS, and US stand for an original replica-exchange simulation, replica-exchange 
umbrella sampling simulation, and conventional umbrella sampling simulation, respec- 
tively. 

b M, Nt, and L are the total numbers of replicas, temperatures, and restraining potentials, 
respectively (see Eqs. fllBD and (H)). In REUS2 and US2 we set M = N T x L for 
simplicity. We remark that this relation is not always required. For instance, the 16 
replicas could have 16 different temperatures with 16 different restraining potentials (i.e., 
M = N T = L = 16). 

di and kg (£ = 1, • • ■ , L) are the strengths and the midpoints of the restraining potentials, 
respectively (see Eq. (|23"D ). 

d The parameter value 0.0 (0.0) means that the restraining potential is null, i.e., Ve = 0. 
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The purpose of the present simulations is to test the effectiveness of the replica- 
exchange umbrella sampling with respect to the conventional umbrella sampling (REUS1 
and REUS2 versus US1 and US2). The original replica-exchange simulation without um- 
brella potentials (REM1) was also made to set a reference standard for comparison. For 
REM1 replica exchange was tried every 20 time steps (or, 10 fs), as in our previous work 
p3| . For REUS simulations, on the other hand, replica exchange was tried every 400 steps 
(or, 200 fs), which is less frequent than in REM1. This is because we wanted to ensure 
sufficient time for system relaxation after A-parameter exchange. 

In REM1 there are 16 replicas with 16 different temperatures listed in Table 1. The 
temperatures are again distributed exponentially. After every 10 fs of parallel MD sim- 
ulations, eight pairs of replicas corresponding to neighboring temperatures were simul- 
taneously exchanged, and the pairing was alternated between the two possible choices 

For umbrella potentials, the 01 to H5 hydrogen-bonding distance, or "end-to-end 
distance," was chosen as the reaction coordinate £ and the harmonic restraining potentials 
of £ in Eq. (^) were imposed. The force constants, ki, and the midpoint positions, d e , 
are listed in Table 1. 

In REUS1 and US1, 14 replicas were simulated with the same set of umbrella potentials 
at T = 300 K. Let us order the umbrella potentials, Ve in Eq. (p2[), in the increasing order 
of the midpoint value de, i.e., the same order that appears in Table 1. We prepared replicas 
so that the potential energy for each replica includes exactly one umbrella potential (here, 
we have M = L = 14). Namely, in Eq. (^2|) for A = A m we set 

A£? = km , (47) 

where 5k,i is Kronecker's delta function, and we have 

E Kt {q®) = E Q {q®) + V m {q®) . (48) 

The difference between REUS1 and US1 is whether replica exchange is performed or not 
during the parallel MD simulations. In REUS1 seven pairs of replicas corresponding to 
"neighboring" umbrella potentials, V m and V m+ i, were simultaneously exchanged after 
every 200 fs of parallel MD simulations, and the pairing was alternated between the 
two possible choices. (Other pairings will have much smaller acceptance ratios of replica 
exchange.) The acceptance criterion for replica exchange is given by Eq. (|T3"|), where 
Eq. ([T7|) now reads (with the fixed inverse temperature f3 = 1/300/cb) 

A = (3 (V m (qW) - V m (q®) - V m+l (qW) + V m+1 (gM)) , (49) 

where replica % and j respectively have umbrella potentials V m and V m +i before the ex- 
change. 

In REUS2 and US2, 16 replicas were simulated at four different temperatures with 
four different restraining potentials (there are L = 4 umbrella potentials at = 4 
temperatures, making the total number of replicas M = Nt x L = 16; see Table 1). We 
can introduce the following labeling for the parameters characterizing the replicas: 

A m = (T m , Am) — ► A/ ; j = (T/, Aj) . , s 

(m=l,---,M) (/=!,■••, N T , J = 1,---,L) 
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The potential energy is given by Eq. (ff^) with the replacement: m — »• J. Let us again 
order the umbrella potentials, Vj, and the temperatures, Tj, in the same order that 
appear in Table 1. The difference between REUS2 and US2 is whether replica exchange 
is performed or not during the MD simulations. In REUS2 we performed the following 
replica-exchange processes alternately after every 200 fs of parallel MD simulations: 

1. Exchange pairs of replicas corresponding to neighboring temperatures, Tj and Tj+i 
(i.e., exchange replicas i and j that respectively correspond to parameters A/ j and 
A/ +1 j). (We refer to this process as T-exchange.) 

2. Exchange pairs of replicas corresponding to "neighboring" umbrella potentials, Vj 
and Vj + i (i.e., exchange replicas % and j that respectively correspond to parameters 
A/ j and A/ j +1 ). (We refer to this process as A-exchange.) 

In each of the above processes, two pairs of replicas were simultaneously exchanged, and 
the pairing was further alternated between the two possibilities. The acceptance criterion 
for these replica exchanges is given by Eq. (EH3), where Eq. ([T7]) now reads 

A = fa - /3 I+1 ) (Eq (gW) + Vj (qW) - E (f) - Vj , (51) 

for T-exchange, and 

A = Pi (Vj (qW) - Vj (f) - V J+1 (gW) + V J+1 , (52) 

for A-exchange. By this procedure, the random walk in the reaction coordinate space as 
well as in temperature space can be realized. 

We now give the details of the results obtained. In order to confirm that our REM 
simulations performed properly, we have to examine the time series of various quantites 
and observe random walks. For instance, the time series of temperature exchange for one 
of the replicas is shown in Figure 5(a). The corresponding time series of the reaction 
coordinate £, the distance between atoms 01 and H5, for the same replica is shown in 
Figure 5(b). 

We see that the conformational sampling along the reaction coordinate is significantly 
enhanced. In the blocked alanine-trimer, the reaction coordinate £ can be classified into 
three regions |3lfl : the helical region (£ < 3 A), the turn region (3 A < £ < 7 A), and the 
extended region (£ > 7 A). Thus, Figure 5(b) implies that helix-coil transitions frequently 
occurred during the replica-exchange simulation, whereas in the conventional canonical 
simulations such a frequent folding and unfolding process cannot be seen. 

After confirming that the present REUS simulations performed properly, we now 
present and compare the physical quantities calculated by these simulations. In Figure 6 
the potentials of mean force (PMF) of the unbiased system along the reaction coordinate £ 
at T = 300 K are shown. The results are from REM1, REUS1, and US1 simulations. For 
these calculations, the WHAM equations of Eqs. ( PH ) and (|2li| ) were solved by iteration 
first, and then Eq. (|27]) was used to obtain the PMF. 

From Figure 6 we see that the PMF curves obtained by REM1 and REUS1 are essen- 
tially identical for low values of £ (£ < 7 A). The two PMF curves start deviating slightly, 
as £ gets larger, and for £ > 9 A the agreement completely deteriorates. The disagreement 
comes from the facts that the average £ at the highest temperature in REM1 (T 16 = 1500 
K) is < £ >Ti 6 ~ 8-0 A and that the original REM with T-exchange only cannot sample 
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Figure 5: Time series of (a) temperature exchange for one of the replicas and (b) the 
reaction coordinate £ for the same replica as in (a) from the replica-exchange umbrella 
sampling simulation (REUS2 in Table 1). 




Figure 6: The PMF of the unbiased system along the reaction coordinate £ at T = 300 K. 
The dotted, solid, and dashed curves were obtained from the original REM (REM1), the 
replica-exchange umbrella sampling (REUS1), and the conventional umbrella sampling 
(US1), respectively. 
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accurately the region where £ is much larger than < £ >t 16 - These two simulations were 
performed under very different conditions: One was run at different temperatures without 
restraining potentials and the other at one temperature with many restraining potentials 
(see Table 1). We thus consider the results to be quite reliable for (£ < 9 A). 

On the other hand, the PMF obtained by US1 is relatively larger than those obtained 
by REM1 and REUS1 in the region of 2 A < £ < 4 A, which corresponds to the struc- 
tural transition state between the a-helical and turn structures. This suggests that US1 
got trapped in states of energy local minima at T = 300 K. In the region of completely 
extended structures (£ > 9 A), the results of REUS1 and US1 are similar but the discrep- 
ancy is again non-negligible. We remark that at T = 300 K the PMF is the lowest for 
£ 2 A, which implies that the a-helical structure is favored at this temperature. 

We next study the temperature dependence of physical quantities obtained from the 
REM1, REUS2, and US2 simulations. In Figure 7(a) we show the PMF again at T = 300 
K. We observe that the PMF curves from REM1 and REUS2 are essentially identical for 
£ < 9 A and that they deviate for £ > 9 A, because the results for REM1 is not reliable 
in this region as noted above. In fact, by comparing Figures 6 and 7(a), we find that the 
PMF obtained from REUS1 and REUS2 are almost in complete agreement at T = 300 
K in the entire range of £ values shown. On the other hand, we observe a discrepancy 
between REUS2 and US2 results. The PMF curve for US2 is significantly less than that 
for REUS 2 in the region 2 A < £ < 8 A. Note that the PMF curves for US1 and US2 are 
completely in disagreement (compare Figures 6 and 7(a)). 
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Figure 7: The PMF of the unbiased system along the reaction coordinate £ at two temper- 
atures, (a) The PMF at T = 300 K. The dotted, solid, and dashed curves were obtained 
from the original REM (REM1), the replica-exchange umbrella sampling (REUS2), and 
the conventional umbrella sampling (US2), respectively, (b) The PMF at T = 500 K. 
The dotted, solid, and dashed curves were obtained from the original REM (REM1), the 
replica-exchange umbrella sampling (REUS2), and the conventional umbrella sampling 
(US2), respectively. 

In Figure 7(b) we show the PMF at T = 500 K, which we obtained from REM1, 
REUS2, and US2 simulations. We again observe that the results from REM1 and REUS2 
are in good agreement for a wide range of £ values. We find that the results from REM1 
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do not significantly deteriorate until £>llAatT = 500 K, whereas it did start deviating 
badly for £ > 9 A at T = 300 K. The PMF curve for US2 deviates strongly from the 
REUS2 results for £ > 6 A and is much larger than that of REUS2 (and REM1) in this 
region. We remark that at T = 500 K the PMF is the lowest for £ « 6 A and low up to 
£ ~ 8 A, which implies that extended structures are favored at this temperature. 

In Figure 8 we show the average values of the reaction coordinate £ of the unbiased 
system as a function of temperature. The results are again from the REM1, REUS2, 
and US2 simulations. The expectation values were calculated from Eq. We find 

that the average reaction coordinate, or the average end-to-end distance, grows as the 
temperature is raised, reflecting the unfolding of the peptide upon increased thermal 
fluctuations. Again we observe an agreement between REM1 and REUS2, whereas the 
results of US2 deviate. 




Figure 8: Average values of the reaction coordinate £ of the unbiased system as a function 
of temperature. The dotted, solid, and dashed curves were obtained from the original 
REM (REM1), the replica-exchange umbrella sampling (REUS2), and the conventional 
umbrella sampling (US2), respectively. Although the highest temperature in REM1 is 
1500 K, only the results for the temperature range between 200 K and 1000 K are shown for 
REM1. Since the lowest and highest temperatures in REUS2 and US2 are respectively 250 
K and 500 K, only the results between these temperatures are shown for these simulations. 



We now present the results of another example of the multidimensional replica- exchange 
method. This time we consider NPT ensemble of argon fluids, and exchange not only the 
temperature but also the pressure values of pairs of replicas during a MC simulation |f56 



Namely, suppose we have M replicas with M different values of temperature and pressure 
(T m ,P m ). The state x of replica i is characterized by the scaled coordinates and the 
volume and its weight is given by 

W ( x ) = e -/^(J5($M)+P m vM)+tfhiVM . (53) 



The transition probability of replica exchange is then given by Eq. fllcf) , where we now 
have 

A = {J3 m - (3 n ) (E (qM) - E (gH)) + (p m P m - (3 n P n ) (v^ - V«) . (54) 

We prepared M = 64 replicas with Nt = 8 temperature and Np = 8 pressure values 
(M = Nt x Np). We alternately exchanged four pairs of temperatures and four pairs 
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Figure 9: Regions where the acceptance ratios of replica exchange become low (shaded 
regions) for (a) temperature exchange and for (b) pressure exchange in the multidimen- 
sional replica-exchange Monte Carlo simulation of argon fluids. The crosses in the grid 
indicate the values of the set (T m ,P m ) in reduced units. 



The shaded regions in Figure 9 are where the acceptance ratios of replica exchange 
become low (< 20 %). These regions are those where the replica-exchange method fails 
due to the existence of first-order phase transitions. The results of Figure 9 suggest that 
the multidimensional REM enables the simulation to connect regions which cannot be 
reached by one-dimensional REM simulations with only T-exchange or P-exchange. 

We now present the results of MD simulations based on replica- exchange multicanon- 
ical algorithm (REMUCA) j|2| . The Met-enkephalin in gas phase was first studied again. 
The potential energy is, however, that of AMBER in Ref. H[| instead of Ref. [ |48|| . In Ta- 



ble 2 we summarize the parameters of the simulations that were performed. As discussed 
in the previous section, REMUCA consists of two simulations: a short REM simulation 
(from which the density of states of the system, or the multicanonical weight factor, is 
determined) and a subsequent production run of MUCA simulation. The former simula- 
tion is referred to as REM1 and the latter as MUCAl in Table 2. Finally, a production 
run of the original REM simulation was also performed for comparison and it is referred 
to as REM2 in Table 2. 

After the simulation of REM1 is finished, we obtained the density of states n(E) by 
the weighted histogram analysis method of Eqs. ([IT]) and (pKp. The density of states will 
give the average values of the potential energy from Eq. (^), and we found 

= < E >t\= —30 kcal/mol , , , 

= < E > Tm = 195 kcal/mol . ^ ' 
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Table 2: Summary of Parameters in REM and REMUCA Simulations 



Run 


No. of Replicas, M 


Temperature, T m (K) (m = 1, • • • , M) 


MD Steps 


REM1 


10 


200, 239, 286, 342, 409, 


2 x 10 5 






489, 585, 700, 836, 1000 




REM2 


10 


200, 239, 286, 342, 409, 


1 x 10 6 






489, 585, 700, 836, 1000 




MUCA1 


1 


1000 


1 x 10 7 



Then our estimate of the density of states is reliable in the range E\ < E < Em- The 
multicanonical potential energy £W (E) was thus determined for the three energy regions 
(E < Ei, Ei < E < E M , and E > E M ) from Eq. ([45]) . Here, we have set the arbitrary 
reference temperature to be T = 1000 K. 

After determining the multicanonical weight factor, we carried out a multicanonical 
MD simulation of 1 x 10 7 steps (or 5 ns) for data collection (MUCAl in Table 2). In Figure 
10 the probability distribution of potential energy obtained by MUCAl is plotted. It can 
be seen that a good flat distribution is obtained in the energy region E\ < E < Em- In 
Figure 10 the canonical probability distributions that were obtained by the reweighting 
techniques at T — Ti = 200 K and T = Tm = 1000 K are also shown. Comparing these 
curves with those of MUCAl in the energy regions E < E\ and E > Em in Figure 10, 
we confirm our claim in the previous section that MUCAl gives canonical distributions 
at T — Ti for E < E\ and at T = Tm for E > Em, whereas it gives a multicanonical 
distribution for E\< E < Em- 
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Figure 10: Probability distribution of potential energy of Met-enkephalin in gas phase that 
was obtained from the replica-exchange multicanonical simulation (MUCAl in Table 2). 
The dotted curves are the probability distributions of the reweighted canonical ensemble 
at T = 200 K (left) and 1000 K (right). 

In the previous works of multicanonical simulations of Met-enkephalin in gas phase 
(see, for instance, Refs. |57|), at least several iterations of trial simulations were re- 



quired for the multicanonical weight determination. We emphasize that in the present 
case of REMUCA (REM1), only one simulation was necessary to determine the optimal 
multicanonical weight factor that can cover the energy region corresponding to tempera- 
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tures between 200 K and 1000 K. 

To check the validity of the canonical-ensemble expectation values calculated by the 
new algorithms, we compare the average potential energy as a function of temperature 
in Figure 11. In REM2 we used the multiple-histogram techniques (or WHAM) |7|, ||, 
whereas the single-histogram method || was used in MUCAl. We can see a perfect 
coincidence of these quantities between REM2 and MUCAl in Figure 11. 
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Figure 11: The average potential energy of Met-enkephalin in gas phase as a function 
of temperature. The solid and dotted curves are obtained from REM2 and MUCAl, 
respectively (see Table 2 for the parameters of the simulations). 

We have so far presented the results of generalized-ensemble simulations of peptides in 
gas phase. However, peptides and proteins are usually in aqueous solution. We therefore 
want to incorporate rigorous solvation effects in our simulations in order to compare with 
experiments. Met-enkephalin was thus studied by both REM and REMUCA simulations 
in aqueous solution based on TIP3P water model |58| . The AMBER force field of Ref . [49[ 



was used. The number of water molecules was 526 and they were placed in a sphere of 
radius of 16 A. Thirty-six replicas that correspond to temperatures ranging from 200 K 
to 700 K were used. 

The time series of the total potential energy for one of the replicas is shown in Figure 
12. We do observe a random walk in potential energy space, which covers an energy range 
of as much as 2,500 kcal/mol. 

For the REMUCA simulation, the multicanonical potential energy and its derivative 
were obtained by the weighted histogram analysis method from the results of a short 
REM simulation (of 100 psec). In Figure 13 the probability distribution obtained by the 
multicanonical production run of this REMUCA simulatoin is plotted. It can be seen 
that a good flat distribution is obtained in the wide energy range. 

Finally, in Figure 14 we compare the distributions of a pair of main-chain dihedral 
angles (<fr, ip) of Gly-2 and Phe-4 around T = 300 K between gas-phase and in-solution 
results. While the results in gas phase are well localized and sharp, those in aqueous 
solution are distributed more broadly. This suggests that the energy landscape in gas 
phase is much more rugged than in aqueous solution; water considerably smoothes out 
the landscape. We remark that a similar observation was made earlier in Ref. |5T5 . 
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Figure 12: Time series of the total potential energy of Met-enkephalin in aqueous solution 
obtained for one of the replicas from the replica-exchange MD simulation. Corresponding 
times series in the canonical ensemble at temperatures 200 K and 700 K are also shown. 
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Figure 13: Probability distribution of potential energy of Met-enkephalin in aqueous 
solution that was obtained from the replica-exchange multicanonical simulation. 
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Figure 14: Distributions of a pair of main-chain dihedral angles (<fr, i/j) of Met-enkephalin 
around T = 300 K for: (a) Gly-2 in gas phase, (b) Gly-2 in aqueous solution, (c) Phe-4 
in gas phase, and (d) Phe-4 in aqueous solution. 
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4 CONCLUSIONS 



In this article we reviewed uses of generalized-ensemble algorithms for free-energy calcu- 
lations in protein folding. 

We introduced two new generalized-ensemble algorithms which are generalizations of 
the replica-exchange method (REM) (we remark that REM is also referred to as parallel 
tempering). The first one is the multidimensional replica-exchange method (MREM), 
with which we showed that the replica-exchange method is not limited to tempering (or 
temperature exchange) and that we can also exchange parameters in the potential en- 
ergy. One particular realization of this method is replica-exchange umbrella sampling 
(REUS) where we perform tempering and/or the exchange of parameters that charac- 
terize the umbrella potential. The second method is the replica-exchange multicanonical 
algorithm (REMUCA), in which we combine the merits of REM and multicanonical al- 
gorithm (MUCA). 

With these new methods available, we believe that we now have working simulation 
algorithms which we can use for free-energy calculations in protein folding. 
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